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ABSTRACT 

We study a series of A^-body simulations representing elliptical galaxies with central masses. 
Starting from two different systems with smooth centres, which have initially a tri axial config- 
uration and are in equilibrium, we insert to them central masses of various values. Immediately 
after such an insertion a system presents a high fraction of particles moving in chaotic orbits, 
a fact causing a secular evolution towards a new equilibrium state. The chaotic orbits respon- 
sible for the secular evolution are identified. Their typical Lypaunov exponents are found to 
scale with the central mass as a power law L oc m" with s close to 1 /2. The requirements 
for an effective secular evolution within a Hubble time are examined. These requirements are 
quantified by introducing a quantity called effective chaotic momentum ££ . This quantity is 
found to correlate well with the rate of the systems' secular evolution. In particular, we find 
that when Jz? falls below a threshold value (0.004 in our A^-body units) a system does no 
longer exhibit significant secular evolution. 

Key words: stellar dynamics - methods: A^-body simulations - galaxies: evolution - methods: 
numerical - galaxies: elliptical and lenticular, cD - galaxies: kinematics and dynamics 



1 INTRODUCTION 

Schwarzschild's (1979) pioneering method confirmed the existence 
of self-consistent models of elliptical galaxies with exclusively 
regular orbits (Schwarzschild 1979, 1982; Richstone 1980, 1982, 
1984; Richstone & Tremaine 1984; Levison & Richstone 1987). 
This, together with the study of the perfect ellipsoid model (de 
Zeeuw & Lyndel-Bell 1985; Statler 1987), which is fully inte- 
grable, contributed to the formation of a prevailing aspect during 
the '80s that the elliptical galaxies in equilibrium mainly consist 
of stars moving in regular orbits, while stars in chaotic orbits are 
either few or inexistent. For this reason galactic dynamical studies 
were based for many years on a systematic exploration of the stable 
families of periodic orbits in static galactic potentials. 

This point of view was challenged in the early '90s when 
Schwarzschild (1993) found that chaos plays a significant role in 
the structure of systems with cuspy density profiles (p; m ; D oc r~ 2 ). A 
number of observations also indicated that high values of Central 
Masses (CM) exist in galaxies (Crane et al. 1993; Ferrarese et al. 
1994; Lauer et al. 1995; Kormendy & Richstone 1995; Gebhardt et 
al. 1996; Faber et al. 1997; Kormendy et al. 1997, 1998; van der 
Marel et al. 1997; van der Marel & van den Bosch 1998; Magor- 
rian et al. 1998; Cretton & van den Bosch 1999; Gebhardt et al. 
2000). Gerhard & Binney (1985) had shown how box orbits are 
converted to chaotic orbits in systems with strong CM. Later, many 
authors studied in detail this phenomenon and its consequences for 
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elliptical galaxies (Merritt & Fridman 1996; Merritt & Valuri 1996, 
1999; Fridman & Merritt 1997, Valluri & Merritt 1998; Merritt 
& Quinlan 1998; Siopis 1999, Siopis & Kandrup 2000; Holley- 
Bockelmann et al. 2001, 2002; Poon & Merritt 2001, 2002, 2004; 
Kandrup & Sideris 2002; Kandrup & Siopis 2003; Kalapotharakos, 
Voglis & Contopoulos 2004; Kalapotharakos & Voglis 2005, see 
Merritt 1999, 2006 and Efthymiopoulos, Voglis & Kalapotharakos 
2007 for a review). In particular, using the Schwarzschild's method 
Merritt & Fridman (1996) found self-consistent solutions for triax- 
ial configurations in the case of a weakly cuspy profile p oc r~ x . In 
these solutions the fraction of particles in chaotic orbits raised up 
to =i 46%. On the other hand, in the case of strongly cuspy pro- 
files p oc r~- the optimal solutions found had a 60% chaotic orbits. 
However, the solutions in the latter case were not fully stationary. 

It should be noted that, while a solution of the Schwarzschild 
method should ideally represent a stationary galaxy, in practice an 
yV-body realization of such a solution will show some secular evo- 
lution due to a number of reasons, namely a) the solution's response 
density may not match precisely the imposed density b) some of 
the chaotic orbits of the solution may not be 'fully mixed' i.e. they 
may have not reached a full covering of their invariant measure in 
the phase space within the method's fixed integration time, and c) 
by its nature Schwarzschild's method cannot secure the stability of 
even a perfect solution. Thus, questions related to the stability or 
secular evolution of a self-consistent model of a galaxy are better 
answered by the N-body method, which is free of the limitations 
of Schwarzschild's method. 

Many works by the A'-body method have addressed the ques- 
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Figure 1. a) The phase portrait (projections on (|y|,j>) of the 4D Poincare surface of sections (z = 0, z > 0)) at an energy level =: 0.74x (the central potential 
well) in a typical case of our A'-body simulations (see below), before the insertion of the CM. We see the invariant curves of box orbits, b) The same phase 
portrait at ; = 90Ti, ma after the insertion of the CM. The invariant curves corresponding to box orbits have been destroyed. At their place we have a chaotic 
domain and an island of regular short axis tube orbits. As the system evolves the chaotic domain is reduced gradually and the island of SAT orbits grows, c) 
At t = 3Q0Thmct after the insertion of the CM, the system has reached its final equilibrium state. The SAT island prevails. In each of these panels blue dots 
correspond to regular motion while red stars correspond to chaotic motion of one of the A'-body particles, d) Before the insertion of the CM the particle moves 
in a box orbit, e) For r = - 180r/ mK ., after the insertion of the CM the particle wanders in the chaotic domain, f) After r — ISOTf^, the particle is trapped by 
SAT tori and its orbit is converted to regular, of the SAT type. 



tions of the orbital analysis and of the relation between the self- 
consistency and the orbital content of either stationary or evolv- 
ing systems (Holley-Bockelmann et al. 2001, 2002; Contopou- 
los, Voglis & Kalapotharakos 2002; Voglis, Kalapotharakos & 
Stavropoulos 2002; Kalapotharakos et al. 2004; Kalapotharakos 
& Voglis 2005; Muzzio, Carpintero & Wachlin 2005; Jesseit, 
Naab & Burkert 2005; Muzzio 2006; Voglis, Stavropoulos & 
Kalapotharakos 2006). In these works it has been repeatedly 
demonstrated that the insertion of a CM in a triaxial system at equi- 
librium leads to significant secular evolution, resulting in a new 
equilibrium state which is, in most cases oblate (Merritt & Quinlan 
1998; Holley-Bockelmann et al. 2001, 2002; Kalapotharakos et al. 
2004; Kalapotharakos & Voglis 2005). 

Kalapotharakos et al. (2004), Kalapotharakos & Voglis (2005) 
extended the research of Voglis et al. (2002) for a series of models 
with CMs. These works unravelled very clearly the transformation 
of the phase space during the secular evolution that leads the sys- 
tems to a new equilibrium state. A summary of our previous results, 
needed in the sequel, can be made with the help of Fig. 1. Fig. la 
shows the phase space structure in an example of our simulations 
before the insertion of the CM. The phase portrait is almost en- 
tirely filled by invariant curves corresponding to box orbits. The 
big blue dots correspond to the Poincare consequents of a specific 
particle's orbit (called hereafter particle A). Note, that these dots 
are provided by the self-consistent A'-body run. After the insertion 
of the CM practically all the tori of box orbits are destroyed and 
their place is occupied by a chaotic domain and by smaller islands 
corresponding to the well known 'boxlets' (e.g. Merritt 1999; these 



orbits were called High Order Resonant Tube (HORT) in previous 
works of ours). Now, as the chaotic orbits diffuse in the phase space 
they cause a secular change of the system's form, which becomes 
more spherical and less prolate. This change, on its turn, affects 
the phase space favouring Short Axis Tube (hereafter SAT) orbits. 
Fig. lb shows the phase space structure (black points) at t = 90 half 
mass crossing times (hereafter Ti, mcl ). We observe a large chaotic 
domain (for |v| < 0.2) and a big island of stability (for \y\ > 0.2) 
corresponding to SAT orbits. Note, that in the present paper we al- 
ways consider the half mass radius R/ t being equal to unity (R/, = 1). 
The red stars correspond to the Poincare consequents of the orbit 
of particle A in the time interval t = - 180r;,„, c ,. The main effect 
is the following: As the island of SAT orbits grows, the majority 
of particles, initially in chaotic orbits, are gradually trapped in the 
domain of SAT tori and the orbits are converted to regular, of the 
SAT type. Fig. lc shows this conversion for particle A. The second 
row of Fig. 1 shows also the gradual conversion of the same orbit 
from box to chaotic and then to SAT, as it appears in ordinary space. 
The phase space structure in Fig. lc corresponds to t = 300Ti, ma (a 
Hubble time) and the big blue dots correspond to the consequents 
of the orbit of particle A within the time interval t = 180-3507/,,,,,.,. 
At / = 300r A ,„ r , the system has already reached its final equilibrium 
state. At this snapshot the number of particles in chaotic orbits is 
considerably lower than at / = 0. 

In the systems studied by Kalapotharakos et al. (2004) the de- 
struction of box tori after the insertion of the CM results in the 
fraction of mass in chaotic orbits rising initially up to the level of 
50% - 80% (depending on the initial maximum ellipticity of the 
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Table 1. The fractions of the total number of particles in chaotic orbits in 
the various systems for t = (the moment when the CM is inserted). 
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Table 2. The fractions of the total number of particles in chaotic orbits in 
the various systems for t = l5QTi lmct (half a Hubble time). 



models 


m = 


m = 0.0005 


m = 0.001 


m = 0.005 


m = 0.01 


Q 


32% 


78% 


80% 


73% 


56% 


c 


23% 


47% 


48% 


40% 


30% 



system). However, during the secular evolution of these systems 
the fraction of chaotic orbits decreases significantly and at the new 
equilibrium state this fraction ranges between ^ 10% and ^25%. 
The rate of secular evolution depends on the mass value of the CM 
as well as on the configuration of the original system before the 
insertion. In particular, different mass values lead not only to dif- 
ferent fractions of particles in chaotic orbits, but also to different 
distributions of the particles' Lyapunov exponents. 
In the present paper our goals are: 

a) to study the relation between the mass value of the CM, on 
the one hand, and the corresponding values of the Lyapunov charac- 
teristic exponents, on the other hand, of those chaotic orbits which 
are produced due to the insertion of the CM and are responsible for 
the secular evolution. We can immediately state the result of this 
investigation, which constitutes one principal result of the present 
paper: we find that this relation is a power law with an exponent 
close to 1 12. Work on the theoretical justification of such a scaling 
law is in progress, but a preliminary discussion is given in the final 
section of the paper. 

b) to determine the specific requirements so that a system ex- 
hibits significant secular evolution within a Hubble time. Surpris- 
ingly, it was found that, independently of the details of a system, 
the onset of significant secular evolution takes place when a quan- 
tity termed effective chaotic momentum surpasses a threshold value 
(equal to 0.004 in the present paper's A'-body units). This thresh- 
old is global, i.e. the same for all the studied systems. This indicates 
that the effective secular momentum is a quantity possibly related 
to a more fundamental statistical mechanical description of the sys- 
tems under study. 

Section 2 briefly describes the models. Section 3 discusses the 
relation between the value of the CM and the characteristic values 
of the Lyapunov exponents of its associated chaotic orbits. Section 
4 contains the definition of the effective chaotic momentum, that 
measures the ability of a system to undergo secular evolution within 
a Hubble time. Finally, section 5 summarizes the main conclusions 
of the present study. 



Table 3. The fractions of the total number of particles in chaotic orbits in 
the various systems at the final equilibrium state. 
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Figure 2. The evolution of the percentage of chaotic orbits in the system 
(Q, m = 0.01). Initially, we have a high percentage (= 80%) of chaotic 
orbits that decreases with time. When the system reaches the final oblate 
equilibrium state the percentage of chaotic orbits remains stabilized at a 
level =s 20%. The horizontal dashed line denotes the percentage of chaotic 
orbits in the Q system (before the insertion of the CM). 



2 DESCRIPTION OF THE MODELS 

Our study refers to the insertion of a CM in a variety of equi- 
librium A'-body systems. The latter are described in detail in 
Kalapotharakos et al. (2004). Here we describe briefly the most 
important features of these systems. 

Two families of models are examined in the present paper. 
They are produced from two different original systems representing 
smooth centre elliptical galaxies in equilibrium. The two initial sys- 
tems are both triaxial, but nearby prolate. The first system (called Q 
system) has a maximum ellipticity E max = 7 and triaxiality = 0.9 
while the second system (called C system) has a maximum elliptic- 
ity E max a 3.5 and triaxiality at 0.8. Note that the short, interme- 
diate and long axes of each system in the present study is aligned 
along the x,y,z axes, respectively. Voglis et al. (2002) calculated 
the fractions of particles moving in chaotic orbits in these systems. 
They found = 32% and a 23% for the Q and C system, respectively. 
These fractions, as well as the identities of the particles moving in 
chaotic orbits, remain almost unchanged in time, before the inser- 
tion of the CM. 

We insert, in these systems, CMs with a variety of mass values 
and follow numerically their evolution by the A'-body code (Allen, 
Palmer & Papaloizou 1990). We then calculate the fractions of the 
particles moving in chaotic orbits at various snapshots of this evo- 
lution. The CM is assumed to have the following density profile 
(Allen etal. 1990): 



Pa 



GM cm cr 



(1) 



where a = 0.005mR g , m = 



is the relative value of the CM with 



2nr{r 2 + a 2 ) 2 

Man 
M, 

respect to the total galaxy mass M g , and R g is the radius of galaxy. 
The profile Q} yields a cuspy profile r for r <K a. We examined 
four different values of m, namely (0.0005, 0.001, 0.005, 0.01). 

The fraction of mass in chaotic motion, corresponding to the 
integration of orbits in a 'frozen' potential of each system at the 
snapshot t = (the moment when the CM is inserted), is approx- 
imately 80% in all the Q models (for various masses m), and ap- 
proximately 50% in all the C models. These percentages turn to be 
rather irrelevant to the values of m considered (Table[]}. 
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Figure 3. Left column: The distributions of log(L ; ) for the chaotic orbits of the Q family systems at the snapshot t = l50Ti, mc ,, with mass parameters m 
as indicated in the panels. A pronounced maximum appears for the systems with m + 0, which is shifted towards higher values of log(L ; ) as m increases. 
Note, that the area below each curve denotes the corresponding function of chaotic orbits. Right column: Same as in the left column for the distributions of 
log(L ci y). 
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Figure 4. Same as in Fig. 3 but for the C family systems. 



6 C. Kalapotharakos 




In subsequent (isochronous) snapshots, these fractions are 
seen to vary for different values of m. Table[2]shows the fractions 
of particles detected in chaotic orbits at the time t = l50Ti, mcl (half 
a Hubble time) after the insertion of the CM. The main remark 
with respect to this table is that, the systems with higher values of 
m (0.005, 0.01) exhibit smaller fractions of chaotic orbits than the 
systems with lower values of m (0.0005, 0.001) which retain a frac- 
tion of chaotic orbits almost equal to the fractions at the snapshot 
t = (compare Tables [T] and [2} . 

Now, Table [3] shows the fractions of the chaotic orbits for the 
systems with high values of in (0.005, 0.01), when the systems have 
reached a new equilibrium state. The main observation here is that 
in all cases the fraction of chaotic orbits at the new equilibrium 
state is far smaller than the fraction of the chaotic orbits in the 
same system at t = (compare Tables [T] and [3). In fact, the fi- 
nal fraction is even smaller than the fraction of chaotic orbits in the 
same systems before the insertion of the CM. In Fig. 2 we see, for 
example, the evolution of the fraction of chaotic orbits for the sys- 
tem (Q, m = 0.01). The gradual decrease of this fraction is shown 
clearly up to the moment when the system reaches the new oblate 
equilibrium state (t 220). In the case of systems with low values 
of m the corresponding fractions were not calculated because these 
systems need much longer intervals than a Hubble time to reach 
an equilibrium. For example, after a period of 20 Hubble times the 
m = 0.001 system of the Q family appears to be still far from equi- 
librium (see fig. 4 of Kalapotharakos et al. 2004). The times (in 
T hmct imits) of final states of all the systems that have reached a fi- 
nal equilibrium state is given in Table 3, for comparison. 



3 CM AND LYAPUNOV EXPONENTS 

In order to distinguish the orbits into regular and chaotic we use two 
indices: a) The Specific Finite Lyapunov Characteristic Number or 
Lj and b) the Alignment Index (AI) or Smaller ALignment Index 
(SALI) (for more details see Voglis, Contopoulos & Efthymiopou- 
los 1998; Voglis, Contopoulos & Efthymiopoulos 1999; Skokos 
2001; Voglis et al. 2002). The Lj has as unit the inverse of the ra- 
dial period of each orbit T T j and provides a good answer regarding 
whether an orbit can be characterized as regular or chaotic (up to 
a certain level). It can also provide a comparison of the relative 
degree of chaos of two orbits provided that their LjS have already 
reached a stabilized value within a given integration time. However, 
since the Lj is measured by a different unit for each orbit, it can- 



not be used for a comparison of the orbits as regards the efficiency 
of their chaotic diffusion in phase space. For this reason, the L,-s 
are converted to quantities with a common unit for all the particles. 
The Common Unit Finite Time Lyapunov Characteristic Number 
or simply L r „ ; is defined by 

L a j = Lj^ (2) 

l rj 

in units l/T/,„„-, where Ti, mr is the period of radial oscillations of 
the orbits with energy equal to the value of the potential at a ra- 
dius equal to the half-mass radius R h . One finds T hmr ^ 3T hmcl . This 
Lyapunov exponent L cu - S compares the chaotic orbits with respect to 
their ability to undergo significant chaotic diffusion within a Hub- 
ble time. In summary, the value of Lj gives us a measure of the 
overall complexity and chaoticity of the phase space inside which 
the orbits lie, while the value of L cll j yields the exponential rate of 
deviation, measured in absolute time units for each individual orbit. 

Now, Fig. 3 shows the distributions of the values of log(L,) 
axis (left column) or the values of log(L c .„ J ) (right column) at the 
snapshot t = 150r/„„ a for the particles in chaotic orbits of all the 
systems of the Q family. The distributions are normalized with re- 
spect to the systems' fraction of chaotic orbits, i.e., the area below 
each curve denotes the fraction of chaotic orbits. The following are 
observed: 

a) All the distributions of the values of log(L^) for m > 
(Fig. 3 - left column) exhibit a global maximum at values log(L ; ) = 
log(L M ). For m = (Fig. 3a) there appears no conspicuous max- 
imum. On the other hand, for m > there appears a conspicu- 
ous maximum which is, furthermore, shifted to the right as m in- 
creases. By examining many orbits with LjS in an interval around 
these maxima we found that such orbits correspond to particles that 
were initially moving in box orbits (before the insertion of the CM). 

b) The distributions of the values of log(L n(J ) axis (Fig. 3 
- right column) are quite different from the corresponding distri- 
butions of log(L,). Nevertheless, these distributions also exhibit 
global maxima, at values log(L n , M ) not very different from log(L M ). 

Similar remarks can be made for the systems of the C family 
(Fig. 4). It is clear from Figs. 3 and 4 that the values log(L M ) and 
log(L njM ) (for m t 0) depend on m. The mean value of L cll j, de- 
noted by L CU j, also depends on m. The main result now, is that all 
the three dependencies are power laws functions of m (Fig. 5). In 
Fig. 5, the red stars refer to the systems of the Q family and the blue 
dots to the systems of the C family. The power laws of L M (Fig. 5a), 
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Figure 6. Each column shows the projections of different mass components of the (Q, m = 0.01) system on the three principal planes (x — z, y - z, x — y) 
at t = 1507"/, mc .,. Left column: Projection of the particles in regular orbits. These particles constitute approximately 44% of the total mass and they form a 
strongly flattened (oblate) configuration. The orbits are mainly of SAT and boxlet (called 'HORT' in Kalapotharakos & Voglis 2005, see text) type. Middle 
column: Projection of the particles in chaotic orbits with Ljs within a window located around the maximum log(L M ) of the distribution of Fig. 3e. These 
particles constitute approximately 36% of the total mass. They form a triaxial configuration and they were moving in box orbits before the insertion of the 
CM. After the insertion, they are responsible, through their chaotic diffusion, for the secular evolution exhibited by the system. Right column: Projection of 
the particles in chaotic orbits with Ljs outside the window around the log(L;v/) maximum. This component constitutes approximately 20% of the total mass 
and it forms an isotropic (spherical) configuration. The spherical distribution is due to the nearly complete mixing of the particles in the phase space due to 
chaotic diffusion. Such particles no longer contribute to the systems's overall secular evolution. 
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Figure 7. Left column: The projections of the particles in regular orbits of (Q, m = 0.01) system on the three principal planes at ( = 300Ti, mcl . These particles 
constitute 78% of the total mass and form an oblate configuration. Right column: The same projection of the particles moving in chaotic orbits. These particles 
constitute 22% of the total mass and form an almost spherical configuration. The chaotic diffusion does no longer have any consequences on the form of the 
system. 



(Fig. 5b) and L cuj (Fig. 5c) are 

L M cc m si (3a) 
L al M k m n (3b) 
L CU j <x m S3 (3c) 

where Sl = 0.56 ± 0.02, .s- 2 ^_0.57 ± 0.07 and ,s 3 = 0.45 ± 0.03. 
The quantities L M , L cuM and L CU j are different measures of the mean 
level of the Lyapunov exponents of the orbits, which have become 
chaotic, precisely by the insertion of the CM. In conclusion, the 
Lyapunov characteristic exponents of the chaotic orbits scale ap- 
proximately as m l/2 . This property appears here as a numerical re- 
sult, but it probably has some theoretical explanation related to the 
scattering of the orbits by the CM near the centre (see discussion). 



4 EFFECTIVE CHAOTIC MOMENTUM 

In Fig. 5 we see that the systems with m < 0.0050, have an L alj 
close to the value 10~ 2 . These systems do not exhibit appreciable 
evolution within a Hubble time. On the other hand, the systems with 
m > 0.0050, develop a well detectable secular evolution within a 
Hubble time. In that case the L cuj values are close to 10~' 5 . In a first 
approximation, we can say that this difference between the mean 
Lyapunov exponents of the systems with low or high central mass 
value m can explain the faster evolution rate in the latter case, since 
the diffusion of a chaotic orbit is effective within a Hubble time 
(~ 100r (!mr ) only when its L cuj value is L cuj > 10~ 2 . 

However, a more careful study revealed two basic require- 
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Figure 8. The distribution of log(L,-) for the particles in chaotic orbits of the 
system (Q, m = 0.01) at three different times a) t = 0, b) t = 1507/, mc , c) 
t = 3Q0Ti lma . The mass converted from chaotic to regular orbits during the 
secular evolution comes mainly from the area around the global maximum. 



ments which must be fulfilled in order for a system to undergo 
significant secular evolution within a Hubble time: 

(i) the fraction of mass in chaotic motion with high L CU j values 
(L CU j > 10~ 2 ) must not be very small and 

(ii) this mass must have an anisotropic initial distribution in or- 
dinary space. In fact, whenever the mass appears isotropically dis- 
tributed (i.e. close to a spherical distribution) this is an indication 
that the chaotic mixing has already taken place in the phase space, 
i.e. the particles in chaotic orbits have already acquired an almost 
uniform distribution in the phase space. We should stress that a 
spherical distribution of the particles in chaotic motion is not in- 
compatible with the appearance of chaos itself, since the particles 
in regular orbits provide a significant perturbation of the overall 
galactic potential from a spherical potential. 

In order to identify which particles in chaotic motion are most 
responsible for the secular evolution of the systems, we collected 
the particles with log(L^) values in the window log(L M ) - 0.3 < 
log(Lj) < \og(L M ) + 0.3 in Figs. 3 and 4, i.e., the particles within 
roughly one dispersion interval around the maxima of each distri- 
bution. The majority of the particles inside this window are found 
to be particles moving mainly in box orbits in the correspond- 
ing smooth centre system before the insertion of the CM. Up to 
t = 150Tf, mc , these particles preserve to a large extend their initial, 
strongly non spherical, distribution. For example, Fig. 6 shows the 
distributions of the particles in regular and in chaotic orbits for the 
system (Q, m = 0.01) at the snapshot t = l50Ti, mcl . The left column 
shows the projections of the particles in regular orbits on the three 
principal planes x-z,y-z and x - y as indicated in the panels. Ac- 
cording to Table[2]these particles constitute 44% of the total mass 
of the system. This component forms a nearly oblate configuration. 
The middle column shows similar projections, but for the particles 
lying inside the +0.3 window around the value log(L M ) = 10~ L62 
of Fig. 3e. This mass constitutes approximately 36% of the total 
mass and has a triaxial strongly elongated configuration. On the 
other hand, the right column shows the distribution for the remain- 
ing particles in chaotic orbits (outside the corresponding window). 
This mass, which constitutes approximately 20% of the total mass, 
creates an almost spherical background. 

The source of secular evolution of the system (Q, m = 0.01) is 
the mass (in chaotic motion) of the middle column. During the evo- 
lution, the majority of the particles' orbits become regular, mainly 



spherical 




oblate. 



E 2 

Figure 9. The time evolution of the various populations of the system 
(Q,m = 0.01) on the plane of the ellipticities (E 2 ,E { ) = (1-f. 1-f). where 
a, b, c are the short, intermediate and long axes respectively as calculated 
by the moment of inertia tensor of each population. The black dots corre- 
spond to the system (Q, m = 0.01) as a whole, the white circles correspond 
to the regular component and the rectangles and the triangles correspond 
to the chaotic component inside and outside of the window of the log(L,) 
distribution, respectively. The total area within each symbol represents the 
percentage of the associated population with respect to the total number of 
particles. The small arrow in each curve indicates the time arrow. The five 
snapshots that are plotted are at t = (0, 60, 150, 300, 600)7)„, IC .,. We observe 
that all the different populations exhibit an overall evolution towards a more 
oblate configuration (E\ = 0). However, the regular component is always 
far from a spherical distribution (corresponding to E\ = Et = 0), while the 
latter is more closely approached from the start by the chaotic orbits outside 
the window. 



SAT orbits (Kalapotharakos et al. 2004; Kalapotharakos & Voglis 
2005). The particles remaining in chaotic orbits form finally an al- 
most spherical distribution. This appears in Fig. 7 which shows the 
particles in regular or chaotic orbits in the same system but for the 
snapshot t = 300Ti, mcl . The particles in regular orbits constitute 
approximately 78% of the total mass (left column). At this time 
the system has already settled down to a nearly oblate equilibrium 
state. It is immediately observed that the particles in regular or- 
bits (left column) support the form of an oblate spheroid, while the 
particles in chaotic orbits (right column) form an almost spherical 
distribution. 

Now, fig. 12 in Kalapotharakos et al. (2004) gives already a 
hint that the orbits converted, during the secular evolution, from 
chaotic to regular belong to the domain of the distributions of 
log(Lj) near the maximum. Fig. 8 shows a comparison of the dis- 
tributions of log(L;) for the mass in chaotic orbits of the system 
(Q, m = 0.01) for three different snapshots, namely t = (dotted- 
dashed line), t = l50Ti lmcl (solid line) and / = 300r;„„ c , (dotted 
line). The area below each curve corresponds to the fraction of the 
mass in chaotic orbits at each snapshot. The clear maximum ap- 
pearing at log(L M ) = 10~' 62 of the first distribution (at t = 0) is 
smaller in the second distribution (at / = 150r;„„ cr ) and even smaller 
(but still detectable) in the last distribution (at t = 300T hmcl ) at the 
position log(L M ) si 10"' 5 . The decrease of this maximum implies 
that many chaotic orbits were converted gradually to regular orbits. 



10 C. Kalapotharakos 



(Q, m = 0.0005) at t = l50T hmct ^ ift 

Particles in chaotic orbits 

regular orbits (22%) inside a window(49%) 



Particles in 
chaotic orbits 
outside a window(29%) 




Figure 10. Similar to Fig. 6 but for the (Q, m = 0.005) system. The regular component (left column) constitutes = 22% of the total mass and forms a strongly 
elongated structure. The part of the chaotic component that lies around the maximum of the distribution of Fig. 3b (middle column) constitutes 49% of the 
total mass and has a very similar configuration to that of regular component. The chaotic diffusion of this mass is mainly responsible for the secular evolution. 
The remaining part of the chaotic component (right column) has a much more uniform and isotropic configuration. These orbits are already mixed in their 
available phase space. 



Such conversions drive the secular evolution of the system, and they 
are produced mainly from particles in the area of this maximum. 

Figure 9 shows how does the distribution of particles in dif- 
ferent types of orbits affect the geometric parameters, i.e. ellipticity 
and triaxiality, of the (Q,m = 0.01) system, in the course of the 
latter's secular evolution. This figure is produced by calculating the 
time evolution of the moment-of-inertia tensor of the (Q, m = 0.01) 



system as a whole, and separately for the regular orbits, the chaotic 
orbits within the main window of the \og(Lj) distribution, and the 
chaotic orbits outside this window. The plot shows the ellipticity 
Ei = 1 - b/c (vertical axis) versus E 2 = 1 - ajc (horizontal axis), 
where a,b,c are the short, intermediate and long axes respectively as 
calculated by the moment of inertia tensor of all the particles of the 
system (black circles). The white circles correspond to the same 
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Figure 12. The effective chaotic momentum Jz? versus m. The red stars 
correspond to the Q family systems and the blue dots correspond to the C 
family systems. In the case of systems showing negligible secular evolution 
within a Hubble time one point is plotted for the value of Jzf at the time 
/ = 1507), mcr . In the case of systems showing detectable secular evolution 
within a Hubble two points are plotted corresponding to t = l5QTi lmct and to 
the equilibrium state. The systems showing negligible or no secular evolu- 
tion have S£ values lower than — 0.004. The systems with detectable secular 
evolution (m = 0.005, m = 0.01) have initially high ££ values (££ > 0.004) 
and finally low ££ values (J£ < 0.004). 



quantities for regular orbits, the squares to chaotic orbits within 
the main window of the log(L^) distribution, and the triangles to 
chaotic orbits outside the main window. The total area within a cir- 
cle, square or triangle symbol is proportional to the percentage of 
the associated population with respect to the total number of par- 
ticles. The five points per curve correspond to five different time 
snapshots as indicated in the figure. Clearly, the overall evolution of 
the system is towards an oblate state (corresponding to the straight 
line Ei = 0). All the different populations exhibit an overall evo- 
lution towards a more oblate configuration. We notice, however, 
that the regular orbits are always far from a spherical distribution 
(corresponding to E\ = E 2 = 0), while the latter is more closely 
approached from the start by the chaotic orbits outside the main 
window. This further substantiates the phenomena described so far 
in Figs. 6 and 7. In particular, the chaotic diffusion causes the orbits 
inside the main window to fill the entire phase space available to 
them. In the configuration space, this implies that the chaotic orbits 



of any given constant energy fill progressively a larger and larger 
volume within the equipotential surface corresponding to that en- 
ergy (this is in contrast to the box or tube orbits which only fill 
boxes or tubes in such a volume). Furthermore, in typical galactic 
potentials the equipotential surfaces are rounder than the surfaces 
of equal density. Thus, the chaotic diffusion causes the spatial dis- 
tribution of the orbits within the main window to become more and 
more spherical and this is the main factor driving the secular evo- 
lution. 

In order to demonstrate that the maxima of the distributions 
of (Figs. 3, 4) are important as regards the secular evolution of the 
systems, Fig. 10 shows a plot similar to Fig. 6, but for the system 
(Q, m = 0.0005) at the snapshot / = 1507a,„ c ,. The particles in regu- 
lar orbits shown in the left column constitute approximately 22% of 
the total mass and form a triaxial strongly elongated configuration, 
which consists mainly of SAT orbits and boxlets that have survived 
after the insertion of the CM. The particles in chaotic orbits which 
are inside the window (±0.3) around the value log(L M ) = 10~ 232 
(Fig. 3b) constitute approximately 49% of the total mass (middle 
column of Fig. 10). The spatial distribution of this mass is quite 
similar to that of the mass of the left column. The right column 
shows the projections of the particles in chaotic orbits with values 
of log(Lj) outside the ±0.3 window (approximately 29% of the total 
mass). These particles form a more isotropic distribution compared 
to these of the left and middle columns in the same figure. Because 
of the low values of log(L M ), the rate of chaotic diffusion of the 
particles in the middle column is so slow that a macroscopic sec- 
ular evolution was not detected even for times much longer than a 
Hubble time. 

In conclusion, the rate of secular evolution of the systems de- 
pends mainly on the fraction of particles in chaotic orbits which are 
distributed anisotropically (strongly non spherical) and on the mean 
value of the Lyapunov exponents L cu j of the same orbits, which 
measures the ability for chaotic diffusion. 

The fraction of anisotropically distributed mass in chaotic mo- 
A/V„, 

tion can be derived by the ratio , where AN W is the mass in- 

Ntotal 

side the window around the maximum of the log(L^) distributions 
in Figs. 3 and 4 and Ntotal is the total number particles in the system. 
If, now, the mean value of L CU j for the mass inside this window is 
L,,,, a measure of the efficiency of the chaotic diffusion can be ob- 
tained by the quantity 

A2V W 



Ntotal 



E w 



(4) 



hereafter called effective chaotic momentum. The physical motiva- 
tion behind the characterization of ,5? as a 'momentum' is that the 
quantity L„. represents a 'speed', i.e. the speed of growth of de- 
viations in the tangent space to the phase space of the particles' 
motion, while this speed is multiplied by a mass AN W , i.e., the total 
mass of the particles in chaotic orbits causing substantial chaotic 
diffusion. 

Now, Fig. 1 1 shows the value of log(L„.) versus log(m) for all 
the experiments. A power law fit L„ oc rrf with s = 0.57 ± 0.02 
is found, which is similar to the power laws found in the previ- 
ous section. On the other hand, the relation between J/f and the 
rate of systems' evolution is shown in Fig. 12. The red stars cor- 
respond to the systems of the Q family while the blue dots corre- 
spond to the systems of the C family. For those systems exhibiting 
detectable secular evolution within a Hubble time, i.e. the Q and C 
systems with m = 0.005 or m = 0.01, the figure shows two points 
per system corresponding to measurements of at two different 
snapshots namely 1 = 1507), mc ., and / equal to time when a system 
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5/oo (slow) B/20 (slow) B/22 (s1qw) 



Figure 13. A comparison of values of the most important coefficients (B/oo. B/20, -B/22, I = 0...19) of the multipole expansion of the potential (see eq. 13 of 
Kalapotharakos et al. 2004) of the final equilibrium state of the system (Q m = 0.01) when the introduction of the CM is done slowly (values in the horizontal 
axis) or abruptly (values in the vertical axis). The points in these diagrams lie close to the diagonal, meaning that the final states derived by the two simulations 
are almost the same. 



reaches its final equilibrium state (as indicated in the figure). For 
the systems with low values of the CM, Fig. 12 shows only one 
point corresponding to the value of Jzf at the snapshot t = l50T hmct . 

Figure 12 shows that Jzf is a good measure of the rate 
of secular evolution. For example, at t = l50T hmc , the system 
(Q, m=0.001) has a very high fraction of chaotic orbits (80%), but 
the majority of these orbits have low values of Lyapunov expo- 
nents (Figs. 3c,c'). On the other hand, at the same time the system 
(C, m = 0.01) has a much smaller fraction of chaotic orbits (30%) 
but the majority of them have high values of Lyapunov exponents 
(Figs. 4e,e'). Thus, the Jzf value of the system (C, m = 0.01) is 
higher than the Jzf value of the system (Q, m = 0.001), In accor- 
dance with that, the system (C, m = 0.01) has a much faster rate 
of secular evolution than the system (Q, m = 0.001). In the case of 
systems not showing detectable secular evolution within a Hubble 
time (Q and C systems with m = 0.0005 or m = 0.001) the value 
of Jzf is low. A threshold value appears to be 0.004. That is, in all 
the evolving systems Jzf initially takes values higher than the value 
0.004. As the systems evolve, the value of Jzf decreases. When the 
final equilibrium is established, Jzf falls below the value 0.004. 



5 DISCUSSION AND CONCLUSIONS 

This paper deals with a series of A'-body models representing ellip- 
tical galaxies with central masses. The main addressed question re- 
gards a quantitative characterization of the chaotic diffusion caused 
by the insertion of a central mass in systems containing initially 
many box orbits, as well as of the consequences of such a diffusion 
process in the secular evolution and macroscopic features of the fi- 
nal states of the systems under study. The following is a summary 
of our principal findings: 

1) The insertion of the central mass initially converts the 
majority of regular box orbits to chaotic orbits. The fraction of 
chaotic orbits raises from 20% - 30% to 50% - 80%. These 'new' 
chaotic orbits diffuse in phase space changing the form of a system, 
which gradually becomes more spherical (less prolate). Due to this 
change, the phase space is gradually transformed as well. The area 
corresponding to short axis tube orbits grows while the chaotic do- 
main is reduced. As the time goes on, many particles in chaotic 
orbits are trapped gradually by the tori of SAT orbits and their or- 
bits are converted to regular. When the secular evolution ceases the 



system settles to a new equilibrium state (in most cases oblate), in 
which it has a low fraction of chaotic orbits (< 25%). 

One point to be commented here is that, in our experiments, 
the CM was "turned on" abruptly while in previous works of the 
same subject (e.g. Merritt & Quinlan 1998, Holley-Bockelmann et 
al. 2001) the authors introduced a gradual increase of the CM up 
to its final value. A question is whether an abrupt insertion affects 
significantly the numerical results and in particular the final state 
reached by the systems. The systems expected to be more sensitive 
on the abrupt insertion are those with the faster rate of evolution. 
When the developing time scale of the CM is much smaller than the 
system's evolution time scale the difference between an abrupt and 
an non-abrupt insertion should be negligible. In our study we are 
interested in CMs with developing times being a small fraction of 
the Hubble time. Thus, systems with long evolution times fulfil the 
above criterion. The more sensitive system is the (Q, m = 0.01) the 
evolution of which towards the final equilibrium state takes place in 
less than a Hubble time. In order to check the difference by a grad- 
ual or abrupt increase of the CM, we evolved this system through 
a different simulation in which the developing time of the CM was 
taken equal to \t Hu bbie (following the ansatz suggested by Merritt 
& Quinlan 1998) and compared the results with those of our orig- 
inal simulation. We checked both the final equilibria and the times 
needed for their establishment. One can easily check that the fi- 
nal states reached by the system with the CM introduced slowly or 
abruptly are quite similar, the only essential difference being that 
when the CM is introduced slowly it takes some more delay time 
for the system to reach the final state. A persuasive test to see that 
the final systems are quite similar is to compare the gravitational 
potential functions in the end of the two simulations. This was done 
by checking the degree of concentration of the coefficients of the 
multipole expansion of the potential (see eq. 13 of Kalapotharakos 
et al. 2004) of the two systems (namely the coefficients correspond- 
ing to the same radial and angular quantum numbers) towards the 
diagonal. Fig. 13 shows this concentration, which demonstrates that 
the final states are quite similar in the two cases, at least macroscop- 
ically. 

2) We showed that the distributions of the logarithm of the 
Lyapunov exponents measured by the inverse radial period of each 
orbit show conspicuous maxima consisting of those chaotic or- 
bits that were previously boxes (before the insertion of the central 
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mass). These orbits are responsible for the secular evolution. The 
Lyapunov exponents of these orbits are related to the value of the 
cental mass by a power law relation LCN oc m" with s close to 
1/2. The chaotic orbits around the maximum are not fully mixed 
in the phase space, i.e., they have not yet covered uniformly their 
available phase space, while for the chaotic orbits not lying near 
this maximum evidence is provided that such a mixing process has 
already taken place. For that reason, the chaotic orbits around this 
maximum have an elongated spatial distribution, while the remain- 
ing chaotic orbits have a much more isotropic distribution. 

An interesting question being raised here regards whether it is 
possible to find a physical justification for the power law scaling 
oc m l/2 of the Lyapunov exponents of the orbits. Theoretical work 
on this direction is in progress and it will be published in a sepa- 
rate study. We only mention that an answer to this question can be 
provided in the framework of the study of orbits in simplified 3D 
potentials like 




V Keplej 



with incommensurable frequencies oj\,(jJ2,ojt,. In such potentials 
we found that the central mass turns most box orbits to chaotic, with 
Lyapunov exponents following precisely the same scaling, i.e., oc 
m 1/2 as in the simulations with the galactic potential. Fig. 14 shows 
an example of this scaling for orbits with initial conditions taken 
randomly on a number of equipotential surfaces of the potential 

Our preliminary investigation shows that this scaling law can 
be explained by an analytical estimation of the eigenvalues of the 
unstable periodic orbits formed in the chaotic domain, as a func- 
tion of m accomplished by the use of perturbation techniques im- 
plemented in the very neighborhood of the unstable periodic orbit 
solutions. Namely, the scaling of the Lyapunov exponents with m 
is related to the scaling of the eigenvalues of these orbits on m. 

3) The level of the Lyapunov exponents is an important factor 
for one system but it cannot alone determine the rate of the secular 
evolution. The rate of secular evolution is significantly different for 
each system and depends both on the value of the inserted central 
mass and on the original structure of one system. The value of the 
central mass, on the one hand, determines the level of the Lyapunov 
exponents (through the power law relation), which are responsible 
for the diffusion rate. On the other hand, the system's original con- 
figuration determines the fraction of box orbits that are converted 
to chaotic orbits able to drive the secular evolution. 

4) A quantity called effective chaotic momentum is defined, 
that is well correlated with the rate of secular evolution of each sys- 
tem. The effective chaotic momentum yields the product between 
the fraction of particles in chaotic orbits, the chaotic diffusion of 
which tends to change the system's form, and the mean Lyapunov 
exponent of the same orbits. Neither high values of Lyapunov ex- 
ponents nor high fractions of chaotic orbits alone can secure an 
effective secular evolution within a Hubble time. An effective sec- 
ular evolution needs a proper combination of these two factors. In 
particular, when the effective chaotic momentum of a system falls 
below the threshold value 0.004 (in the iV-body units) the secular 
evolution practically stops. 
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Figure 14. The LCN oc m 1 ' 2 relation in the case of the toy model potential 
(5). For each value of log(m) we have plotted the logarithm of the mean 
value of the Lyapunov exponents log(LCA') of an ensemble of orbits with 
random initial conditions in the phase space section corresponding to vari- 
ous energy levels. 
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